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ABSTRACT 

We propose a new method to linearise cosmological mass density fields using higher 
order Lagrangian perturbation theory (LPT). We demonstrate that a given density 
field can be expressed as the sum of a linear and a nonlinear component which are 
tightly coupled to each other by the tidal field tensor within the LPT framework. 
The linear component corresponds to the initial density field in Eulerian coordinates, 
and its mean relation with the total field can be approximated by a logarithm (giving 
theoretical support to recent attempts to find such component). We also propose to 
use a combination of the linearisation method and the continuity equation to find 
the mapping between Eulerian and Lagrangian coordinates. In addition, we note that 
this method opens the possibility of use directly higher order LPT on nonlinear fields. 
We test our linearization scheme by applying it to the 2; ~ 0.5 density field from an 
A^-body simulation. We find that the linearised version of the full density field can 
be successfully recovered on >5 /i~^Mpc, reducing the skewness and kurtosis of the 
distribution by about one and two orders of magnitude, respectively. This component 
can also be successfully traced back in time, converging towards the initial unevolved 
density field at z ~ 100. We anticipate a number of applications of our results, from 
predicting velocity fields to estimates of the initial conditions of the universe, passing 
by improved constraints on cosmological parameters derived from galaxy clustering 
via reconstruction methods. 

Key words: (cosmology:) large-scale structure of Universe - galaxies: clusters: gen- 
eral - catalogues - galaxies: statistics 



1 INTRODUCTION 

The present-day mass density field contains information 
about the fundamental pillars of modern cosmology. It is 
a mixture and cross-talk between the primordial hierarchy 
of correlation functions of fluctuations, the law of gravity 
and the value of cosmological parameters. Unfortunately, 
disentangling all these ingredients and extracting useful in- 
formation about them is not a trivial task. On large scales 
this is still relatively simple; linear theory applies and differ- 
ent Fourier modes evolve independently from each other. In 
fact, thanks to these features, cosmological parameters are 
almost routinely constrained using large-scale galaxy clus- 
tering (see e .g. [Cole et al 



2006; Blake et al. 1 120071 



2005 



Eisenstein et al.|2005l|Hutsi 



Percival et al.||2010| [Blake et al. 



2011). The description of small scales is much more diffi- 



cult; highly nonlinear processes are in place, gravity couples 
perturbations on different scales and additional complica- 
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tions arise from nonlinear galaxy biasing and redshift space 
distortions. 

Different approaches have been proposed to recover the 
primordial, linear and Gaussian, density field on medium- 
or small-scales - a process usually referred to as "Gaussian- 
isation" or "Linearisation". The majority are based on a 
local rank-ordered mapping, where the n-th largest density 
fluctuation in one fleld causes the n-th largest perturbation 
in the other. Examples of this are; the logarithm of the lo- 
cal density field ( Neyrinck et al.|2009[ |201l| [joachimi et al. 
2011), a local Gaussianisation assuming that the primor- 
dial probability distribution function (PDF) of densities is 
known CWeinb erg|1992 Yu et al.p Oll ), and applying linear 
or closely linear filters (by simple Gaussianisation (Neyrinck 



et al. 2011) or more sophisticated Wiener-filtering with a 
wavelet truncation (Zhang et al. 201l| ). Although these 



methods have shown to accomplish their goals (with differ- 
ent degrees of success), they still lack a rigorous theoretical 
motivation and support. Additionally, gravity is a nonlocal 
process, where the evolution of a density ffuctuation is not 
only determined by its amplitude, but it also depends on 
the surrounding tidal fleld ( ,Rimes k Hamilton|2006|[2005t . 
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Therefore, a rank-preserving mapping can not be correct in 
detail. 

In this work we propose and explore another way to re- 
cover the initial density field. Our method takes advantage 
of the fact that an initial, linear field and its gravitationally 
evolved counterpart are not independent from each other, 
but are related through the tidal field tensor. This gives us 
an extra piece of information to constrain and improve their 
mapping. Furthermore, the tidal tensor can be predicted an- 
alytically, and be fully specified by a linear field, using higher 
order Lagrangian perturbation theory (L PT, [Buchert et al. 
1994) JBouchet et al.|[T995l |Scoccimarro| |1998| [Bernardeau 



et al.|2002| ). Putting these ingredients together allows us to 
uniquely identify the initial density field which, evolved un- 
der 2LPT gravity, would give rise to the final density field we 
aim to linearise. We note that similar approaches have been 
explored in the past but with limited success and different 
scopes (see Gramann|1993 Monaco & Efstathiou 1999). 

With this physically motivated Gaussianisation process 
we can interprete rank-ordering mappings, in particular, the 
widely used logarithmic transformation. Explicitly, we find 
that the mean transformation between the nonlinear and 
LPT linearised density fields can be approximated by a log- 
arithmic function, consistent with the solution of the linear 
version of the continuity equation. In addition, we demon- 
strate that the linearised field is a good estimate of the initial 
field, not at its respective (earlier) time but at the present. 
Thus, one needs to trace this linearised field back in time, 
or, equivalent ly, to find the transformation from Eulerian to 
Lagrangian coordinates, if this field is to be used for con- 
strained simulations and/or improved cosmological param- 
eters constraints. One exception occurs on large scales or 
high redshifts, where Eulerian and Lagrangian coordinates 
coincide. 

We validate these ideas by applying our 2LPT linearisa- 
tion method to a density field extracted from a cosmological 
A/'-body simulation. The PDF of the resulting linearised field 
is closely described by a Gaussian function. In addition, this 
field correlates with the high redshift outputs of the simu- 
lation, much more than the original field. The correlation 
increases further when the linearised field is traced back in 
time. All this on scales even as small as ~ 5/i~^Mpc. 

An useful consequence of a Gaussianisation, is that the 
linearised field in Eulerian coordinates can be used as an 
input for Lagrangian perturbation theory and consistently 
predict the associated velocity, displacement and future den- 
sity fields. A comparison between estimations of the dis- 
placement field in Eulerian coordinates and the linearised 
field, as given by the logarithmic transformation, was pro- 
vided in a recent paper (Falck et al.||2011| . The latter field 
could be specially useful for reconstruction of the large-scale 
density field in general, and of Baryonic Acoustic Oscilla- 



tions (BAO) in particular (see Eisenstein et al. 



2007 



Noh 



et al.|[20Q9) [Mehta et al. 2011). In a companion paper we 



show that an accurate estimation of peculiar velocities can 
also be obtained in this way. In this case, the usage of our 
estimation of the linear field yields to results superior to 
those obtained by using a logarithmic transformation (Ki- 



taura et al. 



2011 ). In subsequent papers we will address other 



applications of the linearised density field. 

The paper is structured as follows. In section ^ we 
present the theoretical basis for our 2LPT and for the loga- 



rithmic linearisation. We also derive the equations governing 
the time-reversal of the linearised field. We discuss a prac- 
tical implementation in ^and discuss its performance once 
applied to a A/'-body simulation in Q Finally, we present 
our conclusions. 



2 THEORY 

In this section we recap Lagrangian perturbation theory and 
show how a gravitationally evolved density field can be ex- 
panded into a linear and a nonlinear component. We also 
show that the widely used lognormal transformation gives 
an estimate of this linear component in Eulerian coordinates. 
Finally, we derive the equations to trace a density field back 
in time, allowing us to express the linear component (or 
linearised field) in Lagrangian coordinates, which then cor- 
responds to the actual initial density field. 



2.1 Lagrangian perturbation theory linearisation 

Let us start by considering the mapping between the comov- 
ing coordinates of a set of test particles at two redshifts x{z) 
and g(^o), with z < zq. In Lagrangian perturbation theory 
this relation is expressed via a displacement field, "^{q) (see 
e. g. [Bernardeau et al^|2002 ): 



(1) 



which defines a unique mapping between q and x (usually 
referred to as Lagrangian and Eulerian coordinates). We 
note that such a description of gravitational clustering starts 
breaking down when shell-crossing begins. If we further as- 
sume that the test particles were initially homogeneously 
distributed, then we can write the following mass conserva- 
tion relation: 



p{x,z)dx = {p{zQ))dq. 



(2) 



The inverse of the Jacobian of the coordinate transfor- 
mation defines the overdensity field, S = p/{p) — 1, : 



l + 6{x{q,z))^J{q,zy 



with 



dx 



(3) 



(4) 



Combining Eqs. ([T]) and ^ we obtain an expression for the 
density in Lagrangian coordinates q using a relation for de- 
terminants of matrices and assuming curl-free velocity fields 
(^ = — VG, for a discussion on this see Kitaura et al.|20lT ): 



6{q,z) = |l + V,-*(g,z) 



(5) 



^ - V, • *(g, ^) + p^'^ [6] (g, z) + p^'^ [6] (g, z) , 

where the subscripts q refer to partial derivatives with re- 
spect to q. Here one should note that we expand the inverse 
of the Jacobian. The importance of this will become clear 
below. The second term in Jacobian expansion is given by; 
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,(2) 



(6) 

where we use the abbreviation (3,ij = d^S/dqidqj . The third 
term is: 



^(3)[0] ^det(e,. 



(7) 



Note that the density field at redshift z expressed in 
Lagrangian coordinates, 6{x{q,z)), is fully determined by 
the displacement field. This field in turn can be calculated 
within 2nd order LPT (2LPT), in particular it is given in 
terms of two potentials: 

^{q,z) = -D{z)V,4,^'\q) + D2{z)V,<i>^^\q), (8) 
and consequently; 



From which we get: 

6{x,z) = |1-V. (15) 
- V. • z) + /i^'^ [6] (x, z) + /i^'^ [6] {x,z). 

We have thus found that the Eulerian and Lagrangian 
descriptions are equivalent (Eqs. |5]and 15) when the Jaco- 
bian is expanded in the Eulerian frame and the inverse of 
the corresponding Jacobian is expanded in Lagrangian co- 
ordinateo We will use this result in either formulation and 
leave therefore the coordinate dependence out. One should 
note that the formulations given by Eqs. [5] and [15] do not 
transform the density fields from one frame to the other. 
This point will be further clarified in our numerical experi- 
ments presented in §4.1| 

Integrating Eq. ( |11[ ) we get an analogous expression for 
the full potential: 



@{q, z) = D{z)4>^^\q) - D2{z)^^''> (q), (9) 

where D is the linear growth factor, and D2 the second order 
growth factor given by D2 = aD^ and a ~ —3/7. The linear 
0^^^ and nonlinear potential 0*^^-* are obtained by solving a 
pair of Poisson equations: Vq(l)^'^\q) — 5^^\q)^ v^fYieie 5^^\q) 
is the linear overdensity, and Vgc/)*^^-* (g) = 5^'^\q). 

The term 5^'^^ (q) includes the effects of tidal forces and 
represents the 'second-order overdensity' which is related 
to the linear overdensity field by the following quadratic 
expression (see e.g. |Bouchet et ar]|1995 ): 

S^'\q) = Y: {<P^I(<l)<f^]('l) - [<t>':l^ (<!)?) ■ (10) 

i>j 

Inserting these relations in Eq. ^ we get the desired de- 
composition of the field 



D{z)<p, = D{z)^^''> + (16) 

with 0^1^ = -£)2(^)0(2)[<^(l)] + ^(2)[e(^(l))] + ^(3)[e(^(l))] 

and following operator notation (J)'"'^^ [(j)] = V^^/x^^-* [0] and 
= V"^/i^^^[(/)] with (j) being some field. 
This equation tells us how an evolved gravitational po- 
tential is fully determined by its associated linear potential. 
Therefore, linearising a field then becomes an inversion prob- 
lem, which in section ^we discuss how to solve. 

2.2 Lognormal linearisation 

Here we investigate the lognormal transformation as a mean 
to get an estimate of the linear component of the density 
field. Let us follow [Coles Jones| ( |l99l| and start with 
the continuity equation describing the matter content in the 
Universe as a fluid: 



Siq,z) 



5^{q,z) 



+ S'''^{q,z), 



(11) 



where 6^{q,z) = D{z)6^'^\q) is the linear component of 
the density field and the rest being the nonlinear part 
5^1^(9,2) = (2)5(2) (q) ^ ^(2)[e](q,^) + ^(3)[e](q,2). 

From now on we will also use the following notation for 
short 6d = S/D{z). 

Note that Eq. [TT] is only a function of the coordinates 
q and the redshift z. The full nonlinear density field S{q, z) 
is expressed in Lagrangian coordinates while it naturally 
should be expressed in Eulerian coordinates. 

Let us therefore derive the analougous expression in the 
Eulerian frame. We consider now the inverse transformation 
with respect to Eq. [l] 

q = x-'^{x). (12) 
Mass conservation leads now to the following relation 



Nusser et al.|1991 ) 



with 



l + 5{q{x,z)) ^i{x,z), 



(13) 
(14) 



which can be expanded 



dp 1 , 1 „ 

+ - • V) p + -pV • 16 = . 

at a a 



(17) 



(18) 



We can write this equation in Lagrangian coordinates 
introducing the total derivative 



^ = | + i(«.V)p. 

dt at a 



(19) 



If we also switch to conformal time adr = dt then we 



find 



1 dp 
p dr 



-V • u 



(20) 



As long as we can follow particles (no shell-crossings) 
we may also write the continuity equation as 



^ One should note that this equivalency is not true when the Ja- 
cobian is expanded in the Lagrangian frame and then the inverse 
of that expansion is taken as it is done in [Monaco Efstathiouj 
(1999). 
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ln(l + 5) 



-I 



dr V • It . 



(21) 



One must be especially careful at this point as the di- 
vergence of the peculiar velocity field in the right hand side 
of the latter equation is in Eulerian coordinates and not in 



Lagrangian coordinates (for such an approach see jMatarrese 
et aL||1992). The expansion of this term is not straightfor- 



ward for this reason. 

According to LPT (see Eq. |3| we have yet another ex- 
pression for the logarithm of the density field, which can be 
Taylor expanded 

ln(l + ^) = ln(J), 

= in(i- v.-* + /i('^[e] + /i('^[e]) , 

^^^ + ^+(^^), (22) 

where the quantity 5^ summarises all the higher order 
terms. The decomposition in Eq. [22] can always be done. 
The important point to be noticed is that 5^ is in general a 
nonlocal and nonlinear function of 5^ . Taking the ensemble 
averag^of the previous equation we find that 



since {5^ 



(5+) 
0, and 



M = (ln(l + S)} , 



(23) 



the overdensity 6 and a scaled peculiar velocity given by 
V = u/D = u/{fHD) 



dD 



dD 



+ V - ((l + ^)i;) = 0, 

+ (l + ^)V-i; + (i;- V)^ = 0. 



(25) 



We could try to directly integrate this equation back 
in time computing in each iteration the peculiar velocity 
field from the updated density field. However, let us de- 
rive a formulation of the continuity equation which can 
be better compared with previous works and ensures time- 



reversibility. Following Gramann ( 1993 ) we define the devi- 



ation from linear theory as 

Sgy S/D-Si^pt/D 6/D+V-vwith6i^PT = -DVv. 
We can then rewrite Eq. [25] as 

^ - ^ - D {{V - vf + - v)+5gy+D Vidgyv) = . 

(26) 

Under the assumption that flows are irrotational: 
{V -vf + {v-V)V -v ^ ^V^v^ + 26^^^ [0,] (with V = 
—V(j)v), Eq. 26 is simplified to 



06 6 
dD D 



D ( ^ V'^;' + 2/i('^ [0,] ) +5g^+D V-{5g^v) = . 

(27) 



Integrating the latter equation we obtain 



;ln(l + ^)-/i. 



(24) 



In this way we have demonstrated that first order Tay- 
lor expansion of the higher order corrections are given by 
the mean field: {5^) — fi. Note, that in reality the term 
5^ = 5^{x^z) will not be a homogeneous field. In order 
to improve this one has to make higher order expansions. 
The importance of computing the mean field /x was espe- 
cially emphasized in Kitaura et al. ( 2010 ). A way to compute 
this field from the linear field 6 was presented in [Kitauraj 



et al. (2012). Note that the lognormal transformation will 



be equal to minus the divergence of the displacement field 
5^ — — V • ^ only in linear Lagrangian perturbation theory. 



2.3 Time- reversal evolution equations: the 
Eulerian-Lagrangian approach 

In this section we investigate different formulations of the 
continuity equation which permit to trace the structures 
back in time. In particular we find an expression which shows 
how the linear component can be iteratively traced back in 
time. For the derivation of such an equation we will combine 
the results from Lagrangian perturbation theory (based on 
the equation of motion) with an Eulerian formulation of the 
continuity equation. 

Let us express the continuity equation as a function of 

^ The ensemble average is taken over all possible linear fields (in 
Eulerian coordinates) (...} = (... Assuming a /azr sample 

the ensemble average is reduced to a volume average in Eulerian 
coordinates (...) = {• • •)x- 



dD 
with 5q 



1 2 



20^'^ [0.] + 



D^ 



, + V"V-(^^.i;) = 0, (28) 



2.3.1 2nd order continuity equation 

Let us consider only terms up to second order, i. e. ne- 
glecting terms involving 0[D^). The velocity is given by 



(1) 



v'"^ = -V0^^^ + Hereafter the numbers 

in brackets denote the order of the expansion. Accord- 



ingly, the gravitational potential is given by 



(1) 



+ 



*^^'*]. Note, that the peculiar velocity u 
^[1] = in the quadratic term in 

the continuity equation. However, it will have a second order 
contribution in the deviation term 



needs to be linear: 



5PI 



f 



D 



D 



^(2)U(1)1 



(29) 



Putting all together we get 



dD 



1 + 



/ 



D2 

L»2 



a(2)u(i)1 



(30) 

If we neglect the contribution of 2LPT {D2 = 0) we 
recover the formula derived by Gramann (1993). Neglecting 



tidal forces (0^^-* [0*^^-*] = 0) we get the formula by 
Dekell (fT992l 



Nusser & 



© 0000 RAS, MNRAS 000, 000-000 



Linearisation with Cosmological PT 5 



2.3.2 Higher order continuity equation 

To go beyond the Zeldovich approximation in the veloc- 
ity term, say to 2nd order in Lagrangian perturbation the- 
ory, one needs to consider at least 4th order terms in 
the continuity equatio n. T he gravitational potential can be 

a[6] = 



written according to ^2.1 



as 



(1) 



D2 
D 



,(2)U(1)1 



+ 



D \^ ^[O] + 0^^^[O]) with the symbol [6] indicating that 

the 6th order is incomplete. Note that the term 0^^-* [0] = 
li^^\<p] includes sixth order terms involving . How- 
ever a proper sixth order formulation would require includ- 
ing third order Lagrangian perturbation theory. To obtain 
the 4th order equation one would need to truncate that term. 
The deviation term is correspondingly given by 



(31) 



Finally, the continuity equation yields 

, 2 



dD 2 



i(.'^l)'-2,^(^)[,^?l] + l4^+V-V.(<5l 



(32) 



We should mention here that more adequate integration 
solvers are possible which are mass conserving (solvers for 
hyperbolic partial differential equations). However, for the 
studies we are performing in this work the simple scheme 
presented above is adequate. One can notice that the form of 
the continuity equation as given by Eq.[34]is time-reversal as 
it remains invariant under the transformation: (j)'^^^ —cj)'^^^ 
and D -D. 



4 RESULTS 

We carry out our numerical experiments using the Millen- 
nium Run. This simulation tracks the nonlinear evolution 
of more than 10 billion particles, in a box of comoving 
side-length 500/i~^Mpc (Springel et al. 2005). In particu- 
lar, we consider the simulation at different redshifts {z — 
0, 0.5, 1, 2, 3, 127) gridded with nearest-grid-point (NGP) on 
a 256^ mesh. To iterative ly solve the combined Eulerian- 



Lagrangian set of equations described above, we have de- 
veloped a parallel code that uses Fast Fourier Transforms to 
evaluate Laplacian operators and a finite differences method 
for divergence operators. We have dubbed this code as Kl- 



3 NUMERICAL SOLUTION SCHEMES 

Here we present our numerical approach to iteratively solve 
Eqs. pT]and[32l Note that the first equation determines the 
linear component of the field: (j)g[x^ z) (j)^^\x^ z) (the ar- 
row indicates that 0^^-* is calculated from (j)g) and the sec- 
ond equation traces that component back in time yielding 
an estimate of the full component to an earlier cosmic time: 
(j)'^^\x^z) (j)g(x^z + Az) (where Az < in our case of 
study and in this case the arrow indicates that 0^ at an 
earlier time is computed from (j)^^^). 

(i) <t>g{x,z)^4>'-^\x,z) 



We propose to solve Eq. (16) iteratively by updating the 



nonlinear component which depends on the linear potential 



(33) 



with Tg being a scale at iteration i which stabilises the solu- 
tion. Here we use a Gaussian filter with decreasing smooth- 
ing radii. 

(ii) (l)^^\x^z) (j)g{x^z + Az) To compute the time- 
reversal solution we follow [Nusser Dekel ( 1992 ); Gramann 
( 1993 ) and integrate the equation with finite time differences 



(34) 



D21 



(2)r 



4.1 Linearisation 

First, we show how Eq.^^can be used to decompose the full 
nonlinear density field into a linear component and a non- 
linear one. Since LPT breaks down when shell-crossing be- 
comes dominant, we have to smooth the density field to sup- 
press the power on small scales. We apply here a Gaussian- 
kernel with different smoothing radii 5 and 10 /i~^Mpc. We 
also note that the operation of convolution does not com- 
mute with the linearisation. Therefore, we need to ensure 
that this does not seriously affect our results by comparing 
with the true linear field, i. e. with the initial conditions 
of the simulation as we show below. In the upper panels 
of Fig. [1] we solve Eq. [iT] forwards given a linear density 
field taking the first snapshot of the simulation at z = 127 
which we define as the linear component in Lagrangian co- 
ordinates: 6^{q,z = 0) = ^2^°^^(cc = 2; = 127) (middle 
panel) and computing the nonlinear component shown in 
the right panel S^^{q, z = 0). Adding both components we 
get an estimate of the full nonlinear density field at 2; = 0: 
S{q, z = 0). We can see that the nonlinear component is pos- 
itive both in the high and the low density regions in such 
a way that the peaks get more clustered as can be seen in 
the left panel. On the contrary, the voids become less deep. 
This effect is only apparent since the linear component has 
been multiplied by the relative growth factor as explained 
above. However, one should note that all the quantities are 
in the same (Lagrangian) coordinates. The panels in the sec- 
ond row of Fig.[l]show analogous plots but starting from the 
full gravitationally evolved overdensity field at z = on the 
left 6{x,z — 0). Here the linear and nonlinear components 
are computed by numerically solving Eq. [TT] as described 
in §3. Both upper and middle sets of panels look very simi- 
lar, however, a careful inspection shows that the structures 



Kinetic GENeration of initial conditions (in Japanese: origin). 
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6(q,z = 0) 



I 



S^{q, z = 0)[= (52^°^^(a3 = q,z = 127)] 

-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 




6(X,Z = 0)[= (5Nbody(3,^^ ^ 0)] 



6^(x,z = 0) 



6^^{x,z = 0) 



I 



6(q,z = 0) 

_^Nbody(3,^^ ^ 0) 



100 200 300 400 500 



S^{q, z = 0)[= = q, z = 127)] 

-6^{x,z = 0) 



100 200 300 400 500 



X [h-^ Mpc] 




S^^{q,z = 0) 
-5^^{x,z = 0) 



100 200 



400 500 



X [h-^ Mpc] 



Figure 1. Slice through the density field of the Millennium Run after Gaussian smoothing with rg = 10 Mpc. Upper left panel: 
forward solution of Eq. 6 taking as the linear field the Millennium Run at z = 127. Upper middle panel: Millennium Run at z = 127. 
Upper right panel: nonlinear component corresponding to the field in the middle panels. Central left panel: Millennium Run at z = 0. 
Central middle panel: iterative solution of the linear component. Central right panel: nonlinear component. Lower panels: differences 
between the corresponding fields in the upper panels and in the central panels. 
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Figure 2. Matter probability distribution function (PDF) and corresponding skewness and kurtosis for 5]^ : S^^^^^{x^z = 0),ln(l + 
^Nbody^g^^ 2; = 0)) — /i, 5^(£C, 2 = 0) and 5^^{x, z = Upper panels show the decomposition into a linear and a nonlinear component with 
an initial smoothing of rg =5 Mpc (left) and 10 Mpc (right), black: total field, red: linear component, blue: nonlinear component, 
green dashed: lognormal transformation. Lower panels show subsequent steps demonstrating the convergence of the linearisation process. 
Corresponding skewness S and kurtosis K are also indicated. 



are shifted. This is more clearly shown in the lower panels in 
which the differences between both corresponding panels are 
shown. The reason for the shift is that while the upper pan- 
els show the different components in Lagrangian coordinates 
the middle panels show them in Eulerian ones. 

The upper panels in Fig. |2] show the decomposition 
of the fields into a linear and a nonlinear component in a 
more quantitative way for two smoothing scales; 5 and 10 
/i~^Mpc. We show the PDF for the matter in the simula- 
tion at 2; = is shown (black line), the corresponding linear 
(red line) and nonlinear (blue line) components calculated 
with LPT (red line) and the lognormal linearisation (green 
line). We can see that the linearised fields are closely Gaus- 
sian distributed with low skewness {S) and kurtosis (iC), 
whereas the full field, and even more dramatically the non- 
linear component, have considerably large values for S and 
K. A careful inspection of the plots shows that the non- 
linear component does not have a symmetric PDF. This is 
better shown in Fig. [3] The lower panels show the conver- 



gent behaviour of our numerical scheme, demonstrating its 
stable approach to the a solution with progresively smaller 
skewness and kurtosis. 

To further see the effects of the LPT and lognormal 
linear mappings we compute the cell-to-cell correlation be- 
tween the simulation at z = 0: 8^^"^^ {x^ ^ = 0) and the lin- 
ear component 5^{x^z — 0). This can be seen in the upper 
left panel of Fig. [3] We find that the relation between both 
fields is highly nonlinear and that the lognormal mapping 
is in good agreement with the LPT linearisation. However, 
in the LPT case we see a scatter showing that the relation 
is nonlocal. We can see the non- Gaussian nature of the full 
nonlinear field in the x-axis, starting with an overdensity 
5^—1 and reaching moderately large overdensities ^ > 6. 
The linearised field, shown in the y-axis, presents overden- 
sities in the range —2<5<2. The comparison between 
the simulation at 2; = 0: d^^"^^ {x^ z — 0) and the simula- 
tion dX z — 127: 5^^^^ {x — q^z — 127) shows a similar 
relation with a larger scatter. This is due to the fact that 
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^Nbody ^x,Z = Q) (^Nbody (3,^ ^ ^ Q) 

Figure 3. Cell-to-cell comparison after Gaussian smoothing with rg = 10 Mpc between the matter field (5^^°^^(a?, z = 0) of the 
simulation at z = and Upper panels: Left: the iterative solution of the linear component at z = in Eulerian coordinates Right: 
the nonlinear component in Eulerian coordinates, Lower panels: Left: the simulation at z = 127 representing the linear component in 
Lagrangian coordinates, Right: the difference between the simulation at z = and the field at z = 127. The green curve represents the 
lognormal transformation. Dark colour code indicates a larger number of cells and light colour code a lower number. 



apart from the gravitational effects described in the upper 
panels of Fig. ^ there is a transformation from Lagrangian 
to Eulerian coordinates. It is remarkable how well the log- 
normal transformation traces the mean mapping between 
both fields. Additionally, the right panels in Fig. |3] show the 
corresponding nonlinear components. Here we can see how 
the nonlinear field gets positive both in the underdense and 
in the overdense regions compensating for the overestima- 
tion of the deepness of voids in linear theory and largely 
increasing the power in the high density regions. 

4.2 Evolving the linear component back in time 

The purpose of this section is to show that the linear com- 
ponent can be translated from Eulerian to Lagrangian co- 
ordinates. We will make such a demonstration by solving 
Eq. [32] as presented in ^ In the numerical experiments of 
this section we take a starting redshift of z = 0.5 which 
compensates for thte high value of ag employed in the MS 
( Angulo White|[20TQt . 



The results for different redshifts are shown in Fig. [4] 
The left panels show the cell-to-cell comparison between the 
simulation at z = 0.5: 6^^^^^ (x, z = 0.5) and the simula- 
tion at different redshifts S^^^"^^ (x, zj) with zj = 1,2,3. 
One can see in these plots how the relation between the 
fields gets increasingly more biased as expected. The central 
panels show the nearly unbiased cell-to-cell correlation be- 
tween the simulation at different redshifts 5^^^^{x, zj) and 
the time-reversal reconstruction of the full nonlinear density 
field at the same redshift 6^^^{x,z — zj). This demon- 
strates the success in recovering the full nonlinear field at 
scales of 10 h~^Mpc. The right panels show the linear com- 
ponent S^i^prj.(^x, zj) and the tight correlation with the ac- 
tual initial conditions from the simulation. This correlation 
becomes larger with increasing smoothing scale. We have 
denoted the time-reversal reconstructed fields with the su- 
perscript ELPT standing for Eulerian-Lagrangian perturba- 
tion theory due to the combination of both approaches. We 
should note at this point that we have tried the less time- 
consuming approach of estimating the linear field at each 
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52^°^y(a^,^ = 0.5) 5^^''^^{x,z = Zj) S^iq,z = 0) = = 127) 




Figure 4. Left panels: cell-to-cell comparison after Gaussian smoothing with rg = 10 Mpc between the simulation z = 0.5: 
^Nbody^^^^ _ Q^^^ ^j^^ simulation at different redshifts S^^^^ {x^ zj). Middle panels: cell-to-cell comparison between the simulation 
at different redshifts 5^^^^ {x^ zj) and the time-reversal reconstruction of the full nonlinear field at the same redshift 6^^^ (x, z = Zj). 
Right panels: cell-to-cell comparison between the simulation at z = 127 S^{q, z = 0) = 6^°^^ (z = 127) and the linear component of the 
reconstruction at different redshifts: ^-^prj. {x, zj). Note that Zj runs for the following redshifts zj = 1, 2,3. 



time-step with the lognormal approximation. Sorrowfully, 
systematic errors propagate yielding significantly poorer so- 
lutions. We analyze this issue in more detail in [Kit aura et al.| 
(20TT1). 



In Fig. [5] we show the performance of various grid 
based methods to recover the initial conditions including 
the Eulerian-Zeldovich approximation on the left (N usser fc| 
jPekel 1992) , Gramann| ( p^ 93) in the middle and the one pre- 
sented in this work on the right (see {2.3 for a derivation of 
the schemes). Here we integrate the ELPT equations back 
in time up to a redshift of z — 10. For higher redshifts we 



do not observe an appreciable shift in the structures, nei- 
ther an improvement in the correlation to the initial field. 
Moreover, getting stable solutions of the linear component 
becomes more difficult, since the field we are trying to lin- 
earise is already quite linear. The bare eye inspection of the 
plots already shows that the structures are less smooth, the 
voids deeper, and the peaks better confined with increasing 
order of the continuity equation (see upper panels). As a 
consequence, the difference between the reconstructed fields 
and the initial conditions becomes smaller (see middle pan- 
els) . From the comparison between the upper and the middle 
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Figure 5. Upper panels: slices through the reconstructed initial conditions after Gaussian smoothing with rg = 10 Mpc using Left: 
Zeldovich approximation, Middle: Gramann approximation, Right: this work. Middle panels: difference fields between the reconstruction 
and the actual initial field. Lower panels: cell-to-cell comparison between the reconstruction and the actual initial field. 
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Figure 6. Upper left panel: slice through the simulation at z = 0.5 with Gaussian smoothing of rg=5 h ^ Mpc. Upper right panel: 
reconstructed initial condition. Lower panels: difference fields between the corresponding fields in the upper panels and the simulation 
at 2 = 127. 



panels we can conclude that the largest differences are in the 
high density regions. The accuracy of the reconstruction is 
assessed in a more quantitative way in the cell-to-cell cor- 
relations (see lower panels). Clearly, the biases present in 
the Eulerian-Zeldovich approach are considerably reduced 
by taking higher order terms in the continuity equation. 

The power of the Eulerian-Lagrangian reconstruction 
of the initial conditions is more clearly shown with smaller 
smoothing scales. Although our approach will break down 
at scales in which shell-crossing becomes important (as it is 
the case of LPT in general) , we find that it is still extremely 



accurate even for scales > 5 /i ^Mpc (this is further demon- 
strated in a companion paper 'Kitaur a et al.|2011t . In Fig.|6] 
we compare a slice through the simulation at z = 0.5 with 
Gaussian smoothing of rs=5 Mpc (upper left panel) 
with the Eulerian-Lagrangian reconstructed initial field (up- 
per right panel). There, we can see how the clustered regions 
move away from each other and the voids become as strong 
as the peaks when going back in time - a signature of Gaus- 
sian fields. The lower panels hints on the correctness of our 
approach by showing the difference between the correspond- 
ing reconstructed fields and the simulation at z = 127. We 
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Figure 7. Cell-to-cell comparison between the simulation after Gaussian smoothing with rg = 5 /i ^ Mpc at 
simulation at z = 0.5, Right: the reconstruction of the initial field. 



127 and Left: the 



can now appreciate how dipoles, caused by the incorrect po- 
sition of large overdensities, are dramatically reduced with 
ELPT. 

In Fig.|7|we show the cell-to-cell correlation for the sim- 
ulation (after Gaussian smoothing with = 5 Mpc) at 
^ = 127 and both the simulation at z — 0.5 (left panel) 
and the reconstruction of the initial field with ELPT (right 
panel). This quantifies the accuracy of our time-reversal re- 
construction. Here the improvement provided by the ELPT 
scheme in capturing the highly nonlinear and nonlocal re- 
lation between the initial and final fields is evident. The 
correlation between the reconstructed and the actual initial 
field gets significantly tighter and closely unbiased. 

We define the true linear field 5^^^ {k) as the one 
given by the first snapshot in the A/'-body simulation 
5'S'°^^k,z,^^): 5"'^(fc) = C°''''(fc.«init) (^init = 127 in 
the case of the Millennium Run). Accordingly, the linear 
power-spectrum is given by: P^^^ {k) = P^^°^^(A;, ^mit)- We 
should note that the recovered initial fields 8^£^{k^z) are 
smoothed. We need thus to deconvolve the fields to compare 
them to the unsmoothed linear field 5^^^ (k). For this reason 

we define a spherically averaged kernel JC{k,z) = 
which permits us to deconvolve the fields: 



pL IN(fc) 
^(fc,2) 



^rec (^, Z) = JC{k, z)SD''{k, z) . 



(35) 



We define a normalised cross-correlation between two fields 
by 



G(fc, 



v/Pi,-(fc,^)v/J^"N(fc) ■ 
(36) 

We compute G{k, z) for the simulation at redshift z = 0.5 
represented by the red curve in Fig. |8] Any reconstruction 
should yield larger values than this curve. We find that the 
lognormal reconstruction (shown by the green line) gains 
information, however, it introduces small systematic devi- 



ations on large scales (deviations from the unity at scales 
/c^ 0.05). A similar result is obtained by linearising the field 
at ^ = 0.5 with LPT (magenta curve) as described in ^ The 
advantage of this linearisation is that it does not introduce 
systematic effects on large scales. We then trace this field 
backward in time with ELPT to different redshifts: z = 1 
(cyan), z = 2 (blue) and z = 10 (black). We can see that 
there is an important gain of information between z — 0.5 
and z = 2. However, this gain becomes rather moderate 
when going to even larger redshifts. We have checked that 
the Baryon Acoustic Oscillations are significantly recovered 
as expected from the cross-correlation results. Here system- 



atic effects due to nonlinear evolution (see Angulo et al. 
2008| ) are corrected by undoing gravitation within 2LPT (the 
original idea was based on the Zel'dovich approximation, see 
Eisenstein et al. 2007). We will present a more detailed study 
based on a large volume A/'-body simulation in a forthcoming 
paper. 



5 CONCLUSIONS 

In this work we have investigated the linearisation of cosmic 
density fields with nonlocal Lagrangian perturbation theory 
and local rank ordering mapping by a lognormal transfor- 
mation. 

Let us summarise the implications of our findings in a 
series of points: 

(i) Linearisation of cosmic density fields generate estima- 
tions of the initial conditions of the Universe in Eulerian 
coordinates, i. e. in the coordinates in which structures are 
located at present. 

(ii) The relation between the density field and its linear 
component is nonlinear and nonlocal. Local mappings like 
the lognormal transformation can introduce fluctuations on 
large scales that not present in the original fields. 

(iii) The linear component in Eulerian coordinates can be 
used to estimate the peculiar velocity field or the displace- 
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Figure 8. Normalised cross-correlation between the simulation 
at redshift 127 and: the simulation at z = 0.5 (red curve), the 
lognormal transformation (green curve) and the ELPT recon- 
structions at z = 0.5,1,2,10 (magenta, cyan, blue and black, 
respectively). 



merit field using Lagrangian perturbation theory. We further 
demonstrate this in a companion paper ( Kitaura et al.|20lT ). 
Note that the use of the lognormal approximation to obtain 
an estimate of the linear displacement field has been inves- 
tigated in an independent recent work ( Falck et al.|2011 ). 

(iv) The linear component is more correlated with the ini- 
tial conditions than the full gravitationally evolved density 
field and has a potential use to better constrain cosmological 
parameters. This has already been pointed out by 'Neyrinckj 
et al. (2011) for the lognormal case 



The LPT linearisation 
should be even more accurate as it takes the nonlocal tidal 
field component into account. This point remains to be fur- 
ther studied and quantified. 

(v) The linear component can be accurately traced back 
in time on large-scales ^ 5 h~^Mpc from Eulerian to La- 
grangian coordinates yielding fields that are more correlated 
with the initial conditions than the Eulerian representation. 
This implies that Eulerian grid-based methods (as opposed 
to particle based methods, seejPeebles 1989; Nusser & Bran^ 

]|2002^ .Eisenstein et al.,, 20071 



chiml [20001 IBranchini et al.|r2002, .Eisenstein et al.,,2007r 
Lavaux et al.|2008| ) could be used to recover Baryon Acoustic 



Oscillations or other physical signals. 

(vi) We have demonstrated that one can compute the 
nonlinear component from the linear component with LPT. 
It was shown in "Kitaura et al. (2012) (see appendix A) how 
to do that in the lognormal approximation even for the case 
in which the mean field is not known. This can be useful for 
various reasons. It is easier to obtain estimates of the linear 
component than of the full nonlinear density field from ob- 
servational data. The reason being that modeling the power- 
spectrum (or two-point correlation function) in the recon- 
struction method is easier than including higher-order cor- 
relation functions (see Kitaura 2010 ) . It was demonstrated 
Kitaura et al. ( 2012| ) how to transform the density field 



into its linear component to apply a Gaussian prior and de- 
termine the power-spectrum iteratively. One could use the 



same concept with more complex relations between the den- 
sity field and its linear component like the one provided by 
LPT discussed in this work. 

In summary, we have shown how to apply higher order 
Lagrangian perturbation theory to gravitationally evolved 
fields and discussed the manifold of applications which can 
be further developed based on this approach. 
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